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Abstract 

Background: Insecticides are an effective and practical tool for reducing malaria transmission but the development 
of resistance to the insecticides can potentially compromise controls efforts. In this study a mathematical model was 
developed to explore the effects on mosquito populations of spatial heterogeneous deployment of insecticides. This 
model was used to identify important parameters in the evolution of insecticide resistance and to examine the 
contribution of new generation long-lasting insecticidal bed nets, that incorporate a chemical synergist on the roof 
panel, in delaying insecticide resistance. 

Methods: A genetic model was developed to predict changes in mosquito fitness and resistance allele frequency. 
Parameters describing insecticide selection, fitness cost and the additional use of synergist were incorporated. 
Uncertainty and sensitivity analysis were performed followed by investigation of the evolution of resistance under 
scenarios of fully effective or ineffective synergists. 

Results: The spread of resistance was most sensitive to selection coefficients, fitness cost and dominance coefficients 
while mean fitness was most affected by baseline fitness levels. Using a synergist delayed the spread of resistance but 
could, in specific circumstances that were thoroughly investigated, actually increase the rate of spread. Different 
spread dynamics were observed, with simulations leading to fixation, loss and most interestingly, equilibrium (without 
explicit overdominance) of the resistance allele. 

Conclusions: This strategy has the potential to delay the spread of resistance but note that in an heterogeneous 
environment it can also lead to the opposite effect, i.e., increasing the rate of spread. This clearly emphasizes that 
selection pressure acting inside the house cannot be treated in isolation but must be placed in context of overall 
insecticide use in an heterogeneous environment. 



Background 

Malaria is one of the most important parasitic infection in 
humans. Several initiatives from the international health 
community in the past decade have lead to an estimated 
drop in malaria associated mortality from around 1 mil- 
lion in 2000 to about 655,000 in 2010 according to the 
WHO [1], although an independent recent study reported 
the decrease to be from 1.82 million in 2004 to 1.24 million 
in 2010 [2]. 



"Correspondence: sbarbosa@liverpool.ac.uk 

Liverpool School of Tropical Medicine, Pembroke Place, Liverpool, UK 



Among the current recommended interventions to 
control the disease is the use of insecticidal nets or 
indoor residual spraying with insecticide to control vec- 
tor mosquitoes populations [1]. A major issue arising from 
the intense deployment of insecticides is the development 
of resistance to the chemical agents [3] . It is a ubiquitous 
problem, regarded as a major hindrance in the control 
of malaria. Furthermore, the use of insecticides is not 
restricted to public health, in fact, around 90% of all insec- 
ticide is deployed in agriculture [4]. This potential spatial 
heterogeneity of insecticide deployment can give rise to a 
mixed environment for mosquito populations. 
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In the past mathematical models have been used to 
inform resistance management practices [5], determin- 
ing the impact of different mosquito control intervention 
strategies including the protection conferred by bed nets 
[6], and, recently, to develop new approaches such as the 
idea of evolution-proof insecticides [7-9]. However, few 
[10] have considered the spread of resistance in a vari- 
able selection pressure context. Consequently a model 
is presented here that considers different niches in an 
environment, that can offer some insights on the impor- 
tance of different parameters and their interactions in 
the dynamics of insecticide resistance. This approach is 
particularly suitable for investigating the impact of a spe- 
cific long-lasting insecticidal net that is being developed. 
Vestergaard Frandsen has submitted a bed net proto- 
type (PermaNet 3.0) for formal evaluation to the WHO 
Pesticide Evaluation Scheme that incorporates insecticide 
plus synergist. 

Synergists are natural or synthetic chemicals, which 
increase the lethality and effectiveness of currently avail- 
able insecticides, but that are nontoxic to insects on their 
own. They block the metabolic systems that would other- 
wise break down insecticide molecules, helping to restore 
chemical susceptibility that would require higher levels of 
the insecticide [11]. For this reason they are proposed for 
use in overcoming metabolic resistance and also to delay 
the manifestation and/or spread of resistance [12]. 

In this bed net, synergist (Piperonyl butoxide - PBO) 
together with the pyrethroid deltamethrin are incorpo- 
rated into the fibres on the roof panel of the net, while 
incorporating only deltamethrin on a lower dosage in the 
side panels. The rationale behind this approach approx- 
imates the "two-in-one" concept for bed nets. Treating 
different parts of bed nets with different insecticides (e.g. 
combining pyrethroid insecticide, applied to the side pan- 
els of the bed net, together with carbamate insecticide on 
the roof [13]) has been suggested to confer advantages 
over the use of insecticides alone [14]. The assumption is 
that foraging female mosquitoes explore an occupied bed 
net from the top downwards (as the warm air and carbon 
dioxide that emanate from the sleeper move upwards ), 
i.e., will land on the roof first and make their way down 
the side panels. Restricting the synergist to the roof also 



allows the sides of the net to be made of a softer and more 
comfortable fiber for the user [15]. 

In this paper a general and flexible model was developed 
by expanding the usual genetic models to account for spa- 
tial and sexual heterogeneities in insecticides exposures. 
Statistical tools as partial rank correlation coefficients, 
logistic regression and classification trees were used to 
explore specific situations of synergist application and to 
uncover the dynamics of resistance. 

Methods 

Model 

A population genetic model was designed that predicts 
changes in mean fitness and resistant allele frequency as 
outcome variables to explore the relative contribution of 
each different environmental niche to the dynamics of 
the population insecticide resistance status. Mean fitness 
assesses the potential effectiveness of control strategies 
at decreasing the population while change of allele fre- 
quency between generations quantifies selection pressure 
for resistance. 

The model is deterministic, i.e., based on the approx- 
imation of an infinitely large population size so that 
stochastic fluctuations of allele frequencies can be 
neglected. Investigations of the changes in allele frequency 
caused by natural selection are based upon the assump- 
tion that selection operates through differential survival 
of the zygote from birth to maturity. It assumes that ran- 
dom mating occurs among all adults pooled across all 
niches, and that progeny are then randomly distributed 
among the niches. Resistance is determined by one allele 
at one locus [16-18] (S: insecticide susceptible allele; R: 
insecticide resistance allele). 

Table 1 defines the fitness of each genotype for each dif- 
ferent niches. It also defines the proportions exposed to 
each niche, that sum to 1, which implies that a mosquito 
can only encounter a single niche in a generation. 

Four niches were considered: 

1- Insecticide free (n): it can be an area either inside or 
outside a household; 

2- Non public-health related insecticide deployment 
(o): typically insecticide use in agriculture and 



Table 1 Model structure: niches, exposure and genotype fitnesses within each niche 



Niches 




Insecticide free 


Non public 


ITN 


ITN + Synergist 


Exposure males 
Exposure females 
Fitness SS 
Fitr^ess RS 
Fitness RR 


1 - (Oimo + 0(mi) 
1 - (Oifo + Oifl) 
1 

1 -hnZ 
1 -z 


Oifo 
1 -(Po 
(1 - (Po) + hoSo 
(1 - (Po) + So 


0(miO - ^m) 

1 - (Pi 
(1 - (Pi) + hiSi 
(1 - (Pi) + Si 


(^mi^m 
(1 - (pi)k 

[(1 -(pi)+hiSi]k 
[(1 -(Pi)+Si]k 
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households. These are deployed outwith of public 
health mosquito control campaigns, and generally 
out of the control of public health officials; The 
subscript 'o' is used for brevity, noting that casual use 
inside the house, e.g. mosquito coils, would also be 
included in this class; 

3- Insecticide-treated bed nets (ITN); 

4- Insecticide-treated bed nets with synergist on the top 
of the net (ITN + Synergist); 

There is likely to be differential exposure to insecticide 
and hence different selection pressure in the sexes, since 
only females feed on humans and are, therefore, the ones 
most likely to enter human habitations and encounter 
insecticides. Consequently, the proportion of mosquitoes 
{a) that encounter each niche was differentiated by sexes 
and genotype fitness was calculated separately. Fitness of 
the genotype SS in the insecticide free niche was con- 
sidered as the reference fitness level, and other genotypic 
fitnesses were measured relative to this fitness, which was 
taken to be 1. 

Different selection (s) and dominance {h) parameters 
were defined for each niche, except for the insecticide 
free. There is by definition no insecticide exposure in the 
insecticide free niche, so s is replaced by z (the cost of 
carrying a resistance allele). Dominance is not an intrin- 
sic property of the alleles, it depends on the environment 
in which they are expressed, thus the differences in dom- 
inance coefficients between niches. High levels of insec- 
ticide may render the resistance allele recessive, because 
only homozygotes survive, while low levels may allow sur- 
vival of both heterozygotes and resistant homozygotes 
rendering the allele dominant [19-21]. 

This is a highly flexible genetic model, that includes 
a baseline fitness level cp for niches where insecticide is 
deployed, that captures the variable effects on fitness of 
being fully susceptible to insecticides [22]. For example, 
setting (fo = (pi = l means SS genotype always killed when 
meeting insecticides while setting (po = 0.9 means 10% 
of SS will survive exposure in the non public niche. It 
also allows the fitness of a resistance homozygote meet- 
ing a ITN to be less than 1 and therefore smaller than a 
susceptible homozygote in a insecticide free niche, reflect- 
ing the fact that a fully resistance genotype may not be 
completely impervious to the insecticide. For example 
setting (po = 0.9, ho = 0.2, 5^ = 0.6 means that 10% of SS 
genotypes, 22% of RS and 70% of RR survive exposure 
in the non public-health related insecticide deployment 
niche. 

Two parameters were included that relate to the addi- 
tional use of synergist: k that quantifies synergist effi- 
ciency and the proportion of mosquitos meeting both 
insecticide and synergist in the bed net. It is assumed that 
synergist exposure is equally efficient across genotypes. 



For example, if the probability of surviving bed net con- 
tact for 55, RS and RR genotypes is 10%, 22%, 70% (see 
above) and /: = 0.1, the proportion surviving bed net plus 
synergist falls to 1%, 22% and 7%, respectively. It would be 
straightforward to include separate k values for each geno- 
type if the synergist impact differed between genotypes. 
Description of all parameters on Table 2. 

Based on the model described on Table 1 the fitness 
{W with appropriate subscripts) across all niches of each 
genotype will be [23]: 

Males, 

^m,ss = 1 - (Olmo + Olmi) + Olmo (1 " (Po) 

+ ami (1 - Pm) (1 - (Pi) (1) 
+ Oimi Pm (1 - (pi) k 

Wm,rs =[ 1 - (olmo + «m/)] (1 hn z) 

+ amo[iX-(Po)+ho So\ 

(2) 

+ ami (1 - Pm) [ (1 - (pi) + hi Si\ 
+ (^mi Pm [ (1 - (pi) + hi Si ] k 

Table 2 Parameters, symbols and subscripts used In the 



construction of the model 

Parameter Symbol 

Dominance coefficient in each niche h 

Fitness cost of carrying a resistance allele z 

Selection coefficient in each niche s 
Baseline fitness level of susceptible homozygote 

in niches were insecticide is deployed (p 
Proportion of mosquitoes encountering a 

particular niche a 
Proportion of mosquitoes encountering ITN 

that also encounter the synergist ^ 

Impact of synergist k 



k=0; synergist completely effective: all mosquitoes 
encountering insecticide plus synergist die 

/c=1 ; synergist completely ineffective: mosquitoes 
encountering insecticide plus synergist die 

at the same rate as those encountering insecticide alone 

(the model tracks survival in different niches 
-Table 1, so the impact is on a reverse scale) 



Subscripts 

Male m 

Female f 

Insecticide free n 

Deployment of insecticide outside the house o 

Deployment of insecticide inside the house / 
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V^m,rr =[ 1 - (dmo + Cimi)] (1 - z) 
+ Qf^o[(l - (po)+So] 
+ ami (1 - fim) [ (1 - (pi) + Si] 
+ (>^mi i^m [ (1 - (pi) +Si] k 

Females, 

^f,ss = 1 - iplfo + + Otfo (1 - ^o) 

=[ 1 - {afo + a^)] (1 - hn z) 
+ afo[{l -(po) + K So] 
+ afi{l-^f) [{I - (Pi) + hi Si] 
+ afiPf [ (l-(pi) + hiSi] k 



(3) 



(4) 



(5) 



^,rr =[ 1 - {(^fo + (^fi)] (1 - ^) + (^fo[iX - (Po) + So] 
+ afi (l-Pf) [(l-(pi)+Si] 

+ afify [ (l-(pi)+Si] k 

(6) 

If resistance allele frequency is pm Ipf and frequency of 
susceptible allele is Iqj, after selection the genotypic 
frequencies will be: 



RRm = 
RSjfi = 
SSffi = 
RRf = 
RSf = 



Wm,rr Pm Pf 

Wm,rs (Pm qf + Pf ^m) 

^m,ss 

Wm 

^f,rr Pm Pf 

Wf 

^f,rs iPm qf + Pf qm) 
Wf 



(7) 



^f,ss qm qf 



Where W are the mean fitness, given as the sum of the 
numerators: 



,rr Pm Pf + Wm,rs (Pm qf + Pf qm) 

,55 qm qf 



(8) 



Wf = Wf^rrPmPf+Wf,rs (Pm qf + Pf qm) + Wf,ss qm qf 

(9) 



The frequency of the resistance allele in males after 
selection, i.e., in the mating pool for the next generation 
{t + 1), is 



Pm,t+1 



,rr Pm 

Pf + 0.5 Wm,rs (Pm qf + Pf qm) 



(10) 



and the corresponding frequency in females following 
selection is 



Pf,M = 



V^f,rr Pm Pf + 0-5 Wf^rs (Pm qf ^ Pf qm) 
Wf 



(11) 



Under this model, the ratio of change of the gene fre- 
quency per generation is given by 

Pm,t+l 



^Pm = 



Pm,t 

Pf,t^l 
Pf,t 



(12) 



(13) 



All simulations started assuming Hardy- Weinberg equi- 
librium (HWE), but genotypes will move away from HWE 
due to differential selection on the sexes. This reflects 
their different degrees of exposure to different environ- 
ments so that the resistance allele frequency will diverge 
slightly in the breeding individuals of each sex and their 
progeny genotypes will no longer be in HWE. Conse- 
quently, to allow for redistribution of resistance between 
the genotypes the chosen census point was at generations 
10-11, based solely on intuition. 

Parameter values 

The subjective part of the analysis lies in identifying 
plausible values and distributions for the parameters in 
Table 2. 

Initial resistance allele frequency value, /7o=0-001> was 
used in all calculations and was selected to reflect the 
initial stages of insecticide resistance (where most of the 
individual are expected to be heterozygotes). There is lit- 
tle field information available for the parameter values 
appropriate for Anopheles gambiae complex and param- 
eter values vary depending on the species and the local 
environment. The range of values and distributions cho- 
sen were very broad to investigate general properties of 
the system; narrow distributions can be used to investigate 
specific situations. 

The proportion of mosquitoes subject to a particular 
niche {a) were randomly selected from a uniform distribu- 
tion but subjected to the constraint that the sum over all 
niches by sex is 1; values were randomly selected from the 
uniform distribution and then divided by the overall sum. 

The proportion of males that meet the ITN, a^i was 
constrained to always be smaller than the proportion of 
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females aji and smaller than 0.2 (less than 20% of the 
males of the population enter the household and con- 
tact the bed net) to reflect the belief that only a small 
proportion of males enter a household since they do not 
seek to blood feed on humans. The proportion of males 
that is expected to contact the top of the bed net, and be 
exposed to both insecticide and synergist is assumed to 
be very small, so we restricted the maximum value of fij^ 
to 0.2. 

Uncertainty and sensitivity analysis 

Simulations to understand the influence of each param- 
eter on the outcome variables (mean fitness and change 
in resistance allele frequency) were performed using latin 
hypercube sampling (LHS) to generate a data set and 
partial rank correlation coefficients (PRCC) calculated to 
provide a quantitative measure of the impact of each 
parameter [24]. LHS techniques were first developed to 
explore the behavior of complex models in economics, 
engineering, chemistry and physics and have been used in 
models predicting the impact of insecticide-treated nets 
on malaria transmission [25]. 

The analysis was performed using R software [26] and 
implementation of LHS using package Ihs, It does not 
allow for the specification of each variable distribution 
beforehand, so sampling was performed assuming a uni- 
form distribution. Once the sample was generated, the 
uniform sample from a column (variable) could be trans- 
formed to the required distribution (Table 3) by using 
quantile functions (using the qtriangle comand in R). 

Table 3 Parameters range of values used In simulations 

Range of values 



Parameter Minimum Peak Maximum Distribution 







0.001 




Constant 


hn 


0 


0.5 


1 


Triangular 


ho 


0 


0.5 


1 


Triangular 


hi 


0 


0.5 


1 


Triangular 


z 


0 


0.5 


1 


Triangular 


So 


0 


0.5 


1 


Triangular 


Si 


0 


0.5 


1 


Triangular 


(Po 


0 




1 


Uniform 




0 




1 


Uniform 




0 




1 


Uniform 




0 




0.2 


Uniform 


af* 


0 


0.5 


1 


Triangular 




0 


0.5 


1 


Triangular 


Oimi 


0 




0.2 


Uniform 


k 


0 




1 


Uniform 



* Females: all female niches; Males: all male niches expect ITN/ITN+synergist. 



A data set of 3,000 replications was generated, with ran- 
dom parameter sets and the corresponding values of the 
outcome variables using equations 8, 9, 12 and 13. Ten 
replicates of this procedure were performed as suggested 
in [24] to investigate the predictive precision of model 
using LHS as the sampling method. This was achieved 
by analysing each replicate separately and verifying that 
results were consistent across ten replicates. 

Allele frequency ratio under two extreme scenarios of 
synergist effectiveness 

Following uncertainty and sensitivity analysis the evolu- 
tion of the frequency of the resistance allele was investi- 
gated. This was achieved by simulating a scenario with a 
fully effective synergist {k = 0), and the other extreme, a 
scenario where encountering the synergist had absolutely 
no effect {k = 1). The dataset consisted of 3,000 individual 
simulations that were run drawing values from Table 3 for 
the parameters. Each simulation was run twice, once with 
k = 0 and another one with k = 1, The ratio between the 
resistance allele frequency in the population in both sce- 
narios at generation 10 = |^^||Ej) indicates how fully 
effective synergists increase (y > 1) or decrease {y < 1) 
the spread of resistance. 

Results included a counter-intuitive outcome that the 
inclusion of a synergist could lead to an increase in the 
rate of the spread of resistance (i.e. y > 1). Further inves- 
tigation of this result was pursed by performing a logistic 
regression with a binary dependent variable (1 if y > 1 
and 0 if y < 1), therefore quantifying how changes in the 
parameters values affect the odds of getting the unex- 
pected outcome j > 1. In this regression only 14 parame- 
ters out of the 16 could be included. The parameters a are 
codependent, because they must sum to unity, so of^o, afo 
were excluded from the regression since they achieve the 
smaller PRCC values (see later). 

Classification trees 

The model has a substantial number of parameters, 16, so 
the logistic regression becomes inefficient when consid- 
ering all possible interactions between them. An alterna- 
tive approach to logistic regression is classification trees, 
that sub-divide the parameter space into smaller regions, 
where the interactions are more manageable. 

Classification trees are used to predict membership of 
cases in the classes of a categorical dependent variable 
(1 if J > 1 or 0 if y < 1) from their input parameters 
and were implemented using an algorithm that grows a 
binary tree [27]. At each internal node in the tree, a test 
is applied to the input parameters to identify the binary 
distinction which gives the most information about the 
class membership. The process is repeated at each result- 
ing node, continuing the recursion until some stopping 
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criterion is reached where it makes a prediction [27]. 
The threshold of complexity parameter (cp) was one of 
the stopping criteria used here, it ensures that any split 
that does not decrease the overall lack of fit by a fac- 
tor of cp is not attempted; it can be preset or estimated 
using cross-validation. Here it was used cross-validation, 
which is a method for validating a procedure for model 
building, without an independent validation dataset. It 
includes any given random divisions of the data into 90% 
learning and 10% test sets [28]. The optimally sized tree 
was obtained by running 10-fold cross-validations on the 
data and by including another stopping criterion, a min- 
imum of 50 observations in a node in order for a split to 
be attempted. 

Results 

Uncertainty and sensitivity analysis 

LHS was used to generate a dataset for sensitivity anal- 
ysis. The procedure was first replicated 10 times so that 
the model predictive precision could be assessed. The 
standard errors (se) and coefficient of variation (cv) of 
the outcome variables between replications are small (se: 
0.001-0.3; cv: 0-0.005), suggesting that the predictive pre- 
cision of the model does not depend on the LHS generated 
dataset. Statistic evidence (t-test, p-value < 0.05 in all 
replications) indicates a difference between sexes on both 
fitness and rate of change of resistance. Parameter sen- 
sitivity was performed to quantify how a change in an 
input parameter value causes a change in the outcome 
variables. Partial rank correlation coefficients calculated 
between each of the input parameters and the outcome 
variables are shown in Figure 1, in black circles in all pan- 
els. This analysis allows to assess the relative importance 
of the parameters in driving resistance and how it affects 
fitness, especially the magnitude of the correlation with 
the synergist, /c. 

The rate of spread of resistance (Figure 1, A and B)) is 
most sensitive to parameters values of selection and domi- 
nance coefficients and fitness cost (5, /z, z). The correlation 
is negative in niches where insecticide is not employed 
and positive when it is present. In males (A), dominance 
and selection coefficients inside the house (/z/, 5/) have lit- 
tle effect on the ratio of change, presumably because only 
a small fraction is exposed. 

The negative correlation between mean fitness and 
baseline fitness levels penalties {(po and (pi) is the strongest 
of all in both sexes. Male (C) and female (D) mean fitness 
are also sensitive to the parameters a: male PRCC coef- 
ficients are positive with mean fitness in all niches and 
females PRCC coefficients are positive in the insecticide 
free niche and negative in the other two. 

Overall, changes in parameter k do not appear to have 
a big impact in the mean fitness of the population. In 
females the parameter k is positively correlated but small 



in magnitude and Pf (the proportion that meet both the 
synergist and insecticide) shows also only a small negative 
correlation. Both k and Pf show no correlation with mean 
fitness in the male population. 

Allele frequency ratio under two extreme scenarios of 
synergist effectiveness 

Figure 2 shows the rate change of allele frequency com- 
paring the extremes cases of a fully effective {k = 0) and 
of a inefficient {k = 1) synergist. In most cases (90%) y 
is smaller than 1, which is intuitively the most likely out- 
come, i.e., resistance spreads slower in the presence of 
the synergist. Effect of synergist in males and females is 
not strictly comparable but is overall similar. Most impor- 
tantly. Figure 2 shows that the difference between the two 
scenarios is small, most of the ratios are between 0.8 and 1, 
implying that the delay in the spread of resistance caused 
by the synergist is not very large. Nevertheless, what was 
not expected was that in approximately 10% of the cases 
the rate of allele spread can be higher when the synergist 
is fully effective (y > 1), Figure 3 shows the predicted fre- 
quency of the resistance allele under different values of k 
(ranging from 0 to 1) in a scenario which j > 1 to illustrate 
the difference in the spread of resistance. As an example, 
at generation 70 the predicted frequency when the syner- 
gist is inefficient {k = 1) is 0.11 and when is fully effective 
(k = 0) is 0.26. 

Logistic regression 

Results of the logistic regression are presented in Table 4. 
Each one unit of increase in the parameter value in ques- 
tion will increase/decrease (+/— signal of the estimate) 
the log odds of the unexpected event {y > 1), the odds are 
the exponentiated values of the estimates {OR = g^^^^^^^^), 
shown also in Table 4. The parameters fy, and a mi 
are not significative (p-values > 0.05) and appear to have 
no impact on the outcome. An increase in the parame- 
ters So, ho, and (po increases considerably the odds off the 
counter-intuitive result (y > 1) and a increase of ajn, hi, 
Sh <Pu olfiy hyi and z slightly decreases the odds. Parameters 
related to the niche where insecticide is not deployed {n) 
do not have much impact on the counter-intuitive result, 
which appears to be governed mainly by the values of the 
parameters in the niche where insecticide is being applied 
for other reasons outside the house (0) and the niche 
insecticide inside the house {i). 

Additionally, to compare these two niches, a series of 
Wilcoxon signed-rank tests were performed. From the 
results of the tests it was determined that all the parame- 
ters in the niche where insecticide was encountered out- 
side the house (o^mo; (^fay ^o> and ho) were significantly 
higher than the equivalent parameters in the niche inside 
the house in the simulations that led to the unexpected 
outcome. 
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A - Ratio of change males 



B - Ratio of cliange females 



o 
o 

DC 
Q. 



o 
o o 



o Original dataset 
O Constrained dataset 



Q /s /s /\ <r> 



o o 



o o 



8 



8 



o 
o 



I I I I I I I I I I I I I I I I I I I 

hp ho hi z So Si cpo cpi k Pf Pm Pf Pm CXfo ttfi Ofp Otmi Otmr 

C - Mean fitness males 





I I I I 

hp ho hi z s, 



I I I I I I I I I I I I I 

Si qpo cPi k Pf Pm Pf Pm Ctfo Ofi afn amo Otmi Otm 

D - Mean fitness females 




hp ho hi z So Si cpo qpi k pf Pm Pf Pm CXfo Otfi OCfp Otmi Otm 



I I I I I I 

hp ho hi z So Si cpo cpi k Pf Pm Pf Pm Wfo OCfj Ofp ttmo Otmi Wm 



Parameters 

Figure 1 Plot of Partial Rank Correlation Coefficients between each of the parameters in the model and the two outcome variables: mean 
fitness and ratio of change of resistance allele frequency, in both sexes (above zero increasing positive correlation, below zero increasing 
negative correlation. The black circles refer to the coefficients calculated using the original dataset and the red diamonds coefficients were 
calculated using a dataset generated with the constrains: afn < 0.2, Sq > 0.5, cpi < 0.8, s,- < 0.5, derived from the classification trees analysis. The 
parameters symbols in the x-axis are defined in Table 2. 



Classification trees 

The classification tree in Figure 4 is a tree pruned to 
avoid overfitting the data, that minimizes the cross- 
vaUdated error [28]. The parameters actually selected by 
the algorithm (shown to have discriminant value) to con- 
struct the tree shown were: hi, ho, Si, Sq, oLfyi, ajo, (pi 
and (fo. The proportion of observations correctly clas- 
sified at each leaf can be used to represent the likely 
proportion of similarly classified observations of unsam- 
pled data at the field conditions defined by that termi- 
nal node [29]. The proportions of classifications on the 
five terminal nodes that predicted the class 1 ranged 
from 0.059 ((j^) to 0.35 {^^)^ Further simulated 
datasets, with higher number of observations, produced 
slightly different trees, but they agree on the parame- 
ters selected for their construction and have the same 
basic structure. 



The parameters most closely associated with the 
counter intuitive results are consistent between logistic 
regression and classification trees. Logistic regression is 
considered to be a more potent method, [30] but the 
schematic nature of the trees provides a clearer under- 
standing of the interactions between the parameters and 
does offer a set of rules to follow in order to try and 
achieve a particular outcome. 

Constrained datasets 

Objective analysis using logistic regression and classifica- 
tion trees led to a subjective conclusion as to why counter- 
intuitive results iy > 1) occur (see later). To validate the 
results new datasets with constrains on the parameters 
based on the tree decision criteria (e.g.: < 0.2, 5^ > 0.5, 
(Pi < 0.8, Si < 0.5), were generated with the expectation of 
reducing substantially the number of observations where 
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Figure 2 Histogram of the ratio of resistance allele frequency at generation 1 0 in the extremes cases of a fully effective {k = 0) compared 
to an inefficient (A- = 1) synergist: y = ^|^|^^^ ). Only values ofy < 2 shown, which constitute 99.3% of the number of simulations. Values higher 
than 1 (y > 1 ) indicate the counter-intuitive result, i.e., that the synergist presence drives resistance faster. 



y > I. From 3,000 simulations produced with the previous 
constrain none resulted in y > 1. To end the analysis, a 
new dataset was generated with the above constrains to 
examine the impact of the different parameters on mean 
fitness when all observations lead to y < 1, Figure 1 
shows in red the PRCC results calculated with this con- 
strained dataset. The estimates values for the original 
dataset (black circles) and the constrained dataset (red 
diamonds) are similar, which implies that the presence 
of the counter-intuitive scenarios where y > 1 did not 
influence the overall correlation between the parameters 
and mean fitness. There was, however, a difference in the 
PRCC between the parameter 5^ (selection in the environ- 
ment with insecticide outside the household) and the ratio 
of change of resistance allele frequency in the male popu- 
lation, it shows no correlation in this dataset while it did 
in the original dataset. The most plausible explanation is 
that male selective pressure occurred mostly in the niche 
with insecticide employed outside. This selection coeffi- 
cients were reasonable high in the simulations that led to 
the unexpected outcome, so it can be considered normal 
that it showed impact in the original dataset and not in the 
new dataset. 

Dynamics of spread 

The dynamics of spread was investigated by checking 
allele frequency at generation 1000 and determining 



resistance status (allele fixed if frequency greater than 
0.99, lost if smaller than 0.001 and at equilibrium if change 
of frequency smaller than 0.001 in the last 50 genera- 
tions) in the original dataset, subset with simulations that 
scored y > 1 and constrained dataset. Results are shown 
in Table 5. 

Discussion 

Protection against vector borne diseases predominantly 
depends upon the usage of insecticides. Different strate- 
gies of delivery, in combination or independently, can 
be enforced while trying to minimise the emergence or 
impact of resistance. This study presents a model where 
mosquitos face an heterogeneous environment of four dif- 
ferent niches, considering the use of insecticides outside 
the household and the existence of insecticide free areas 
(refugia). The model also allows to check the effect on 
resistance spread of new generation long-lasting insecti- 
cidal nets, that incorporate a synergist, reported to have 
improved increased efficacy against pyrethroid-resistant 
malaria vectors [15,31]. It would be simple to add more 
niches, essentially adding extra columns to Table 1, which 
is a major benefit of developing such a flexible method- 
ology. For example, the outside niche could have high or 
low levels of insecticide. Under these conditions of differ- 
ent insecticide concentrations the resistance allele may be 
recessive or dominant respectively, with a huge impact on 
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Figure 3 How synergists affect the spread of resistance: Pre- 
dicted resistance allele frequency in females for different values 
of synergist impact {k = 0, 0.2, 0.5, 0.7, 1) in a specific setting 

{hn = 037 y ho = 0,50, hi = 0,07 yZ = 0,87, So = 0.47, s/ = 0.46, 
(Po = O.lh <pi = 033yPrn = 0,0% = 0.82,0/0 = 0.41,0^ = 0.50, 
afn = 0.10, a^o = 0.38, a^/ = 0.34, = 0.28) that scored 
y > 1. 



rate of resistance [20] . It would also be possible to allow 
for mosquitoes to be exposed to more than one niche, by 
multiplying the fitness in each niche. This would however 
increase the complexity of the model and, as presented, 
the model demonstrates how interesting results can be 
derived from simple approaches. 

Here it was only considered the existence of different 
niches with different use of single insecticides, not allow- 
ing for deployment of a second insecticide with a different 
mode of action. This would require a new model assum- 
ing 2 loci, each encoding resistance to a single insecticide. 
These 2 locus models are simple in principle [32] but in 
the present model of multiple niches it would increase the 
number of parameters significantly (total of 9 genotypes 
for each sex and increase the number of potential niches to 
reflect different combinations of insecticides) and would 
not be beneficial in the current exploration of the effect of 
the synergist. 

The calibration of the model with field data proved to 
be problematic, hence the decision to sweep a range of 
values and check the outcomes. The restriction on male 
proportion inside the house, less than 20%, was set based 
on personal communication since the numbers of males 
that enter households is rarely reported in hut trials, con- 
sequently it must be seen as a rough estimate. This work, 
in spite of the simplicity of the model, is illustrative of 
the advantages of modelling; an overall understanding that 



Table 4 Logistic regression: parameter coefficient 
estimates, standard error, p-value and odds ratio 





Estimate 


Std. Error 


p-value 


OR 


Intercept 


7.614 


0.975 


< 0.05 


2.026e+03 


hn 


-3.215 


0.553 


< 0.05 


4.016e-02 


ho 


7.828 


0.685 


< 0.05 


2.511e+03 


hi 


-8.761 


0.696 


< 0.05 


1 .567e-04 


z 


-3.057 


0.515 


< 0.05 


4.703e-02 


So 


8.109 


0.688 


< 0.05 


3.324e+03 


Si 


-7.645 


0.681 


< 0.05 


4.785e-04 


(Po 


3.797 


0.416 


< 0.05 


4.455e+01 


<Pi 


-6.929 


0.544 


< 0.05 


9.794e-04 




0.190 


0.387 


0.6 


1 .209e+00 




0.378 


1.833 


0.8 


1 .460e+00 


Oif, 


-6.584 


0.753 


< 0.05 


1.383e-03 


Oifn 


-22.538 


1.447 


< 0.05 


1.628e-10 


0(mi 


-2.205 


3.016 


0.46 


1.103e-01 




-0.948 


0.472 


0.04 


3.875e-01 



would not be possible by specific calibration and also 
the emergence of non-intuitive outcomes. Mathematical 
models have been used to expose other so non-intuitive 
results such that indoor residual spray (IRS) of insecti- 
cides in conjunction with bed nets can show antagonism, 
arising via interference of their modes of action while 
it is mostly assumed that the two tools have synergistic 
benefits in reducing malaria transmission [33]. 

Experimental studies have been conducted in the field 
to assess the potential of prototypes of bed nets that incor- 
porate a pyrethroid insecticide in the side panels and the 
synergist PBO plus insecticide on the roof. An experi- 
mental hut trial in Tanzania was not able to demonstrate 
an improvement in efficiency (measured mortality, pas- 
sage through holes and feeding rates) when compared to 
standard insecticide impregnated nets against Culex quin- 
quefasciatus [31] and only moderate performance against 
pyrethroid-resistant A. gambiae M taxon was reported in 
another trial in southern Benin [34]. However, Killeen and 
colleagues [35] noted that the manufacturer of the proto- 
type claims only that the net has greater efficacy than its 
predecessor and that their simulations corroborates this 
claim, which the findings of this work also support. 

From the PRCC analysis is straightforward to infer that 
the success of control campaigns depends mostly on the 
proportion of mosquitoes that encounters insecticides 
and on the fitness scaling factors (p, that define survival 
of SS genotypes when meeting insecticides. These last 
parameters were introduced to emphasise possible dif- 
ferences among niches in the impact on fully susceptible 
SS genotypes when facing insecticides. Unsurprisingly the 
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Figure 4 An example of a classification tree of outcomes y > 1 (class 1 ) and y < 1 (class 0) using all the parameters in the model except k. 

The items displayed in tine nodes in tine tree diagram are: tine criterion for mal<ing tine decision {e.g.-.afn > 0.1 626), tine predicted class for tliat terminal 
node (0 or 1 ) and the number of observations correctly classified to the class versus the number misclassified, achieved by cross-validation (e.g.: 
21 62/84). To proceed through the diagram at a given node move to the left branch if the stated condition is true (yes) and to the right if false (no). 



control measures are more effective the higher their val- 
ues. The inclusion of this parameter was crucial because it 
considers the complexity of fitness, and incorporates the 
differential environmental effects of insecticides across 
different genotypes. 

The present model includes two parameters that relate 
to the effects of application of a synergist in combination 
with a insecticide, /c, the reduce survival due the synergist 
and the proportions of mosquitoes, P (males and females), 
that meet both chemicals. The magnitude of the corre- 
lation between males mean fitness and k and firn is very 
small, presumably due to the small proportion of males 
that is exposed to insecticide in an household, but k and 
P only correlates moderately with females mean fitness 
and ratio of change of allele frequency. This indicates that 
the synergist has a small impact in controlling the pop- 
ulation, but even small values of k will help to recover 
the effect of the insecticide, and possibly this is the main 
contribution of the synergist. Nevertheless adding syner- 
gists to bed nets does decrease the rate at which resistance 
spreads in about 90% of scenarios. 

It was not possible to use the model to investigate the 
overall impact of changes in fitness in the mosquito pop- 
ulation dynamics, which would require a more elaborate 



model incorporating demography (e.g. [36]), and more 
specifically to investigate the effects of targeting mainly 
females, decreasing their fitness more than that of males 
which are generally regarded as determining overall pop- 
ulation regulation. The finding that a situation can arise 
in which having a fully effective synergist in place con- 
tributes to intensify the spread of resistance is the most 
interesting result of this work. The setting in which it 
emerges is very specific, it encompasses a strong selective 
pressure for resistance in the niche outside the household 
(mean So = 0.62 while mean 5/ = 0.41; these and the fol- 
lowing figures come from the simulated subset where y > 
1), most mosquitoes are exposed to insecticides (mean 



Table 5 Dynamics of spread of resistance allele frequency 





% Loss 


% Fixation 


% Equilibrium 


Original Dataset 


49 


30 


21 


Original Subset with y > 1 


3 


72 


25 


Constrained dataset 


8 


79 


13 



Percentage of simulations that lead to loss, fixation and equilibrium of 
resistance allele in the original dataset, subset of the original dataset with 
simulations that scored y > 1 and constrained dataset based on the 
classification tree criteria {afn < 0.2, So > 0.5, (p, < 0.8, 5, < 0.5). 
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for females is 0.16 and 0.44 for males) and there are 
different dominance values in the niches were insecticide 
is deployed (mean ho = 0.61 and mean hi = 0.37). 

In these circumstances the insecticide in the bed nets 
will be largely ineffective and the pressure for selection 
is weak. It seems as if this niche is acting as a refugia 
for susceptibles, that will contribute with their suscepti- 
bles genes for the next generation, therefore decreasing 
the resistance allele frequency. If a fully effective syner- 
gist {k = 0) is present, the fitness of all genotypes inside 
the house will be zero (k affects the 3 genotypes equally, 
so all mosquitoes die irrespective of their genotype) and 
the next generation will be mostly composed by progeny 
of survivors from the niche outside the household where 
selection for resistance was high. An hypothesis is that in 
this particular case the synergist is removing the refugia of 
weak selection in the house thereby magnifying the effects 
of selection for resistance outside the house. 

Males and females showed the same patterns of spread, 
49% converged to fixation, 30% to loss and 21% to equilib- 
rium. Equilibrium in single niche models (such as the sin- 
gle equation approach used in [20]) can only be achieved 
if there is heterozygote advantage [37] which was not 
postulated in the model because dominance coefficient 
h lie between 0 and 1 (Table 3). The balancing effects 
of different selection and dominance acting in different 
environments in the same population seems to be able to 
keep the resistance allele frequency at equilibrium. As an 
hypothetical example, a mutation which is dominant for 
insecticide resistance but has a large, recessive effect on 
fitness. As resistance starts to spread, most R alleles are in 
heterozygotes which resist insecticides and do not pay the 
fitness penalty. As R increases in frequency the proportion 
of RR homozygotes increases, so fitness penalties esca- 
late until an equilibrium occurs. In effect, the marginal 
fitnesses (i.e. the average over all niches) generated by 
Table 1 and Equations 1 to 7 result in the heterozygote 
being the most fit in some simulations. As far as it was 
possible to verify it has not been reported in the field, 
possibly because it is not currently regarded as a likely 
occurrence. 

The three dynamics of spread were predicted in subse- 
quent analysis (Table 5) by estimating allele frequency for 
1,000 generations. In the simulations that scored 3/ > 1, 
3% eventually lead to loss, 72% to fixation and 25% to equi- 
librium. It is overall a worse picture than with the original 
dataset. The choice of parameter values is important, for 
example constraining the dataset reduces the possibili- 
ties of reaching equilibrium, i.e., fixation of resistance was 
much more likely. Analysing only the simulations that lead 
to fixation in the original dataset generated results very 
similar (not shown) to Figure 1 and Table 4. 

These results emphasise a very important fact often 
overlooked in modelling resistance: that it is highly 



dangerous to consider selection in only a single niche, iso- 
lated from other selection pressures, and then extrapolate 
the results from the single niche to the whole popula- 
tion. In this case it seems reasonable to conclude that 
adding effective synergists will reduce selection for resis- 
tance in the household niche because all three genotypes 
are killed. The impact that a fully effective synergist will 
have in disease transmission is a fundamental question, 
that cannot be directly answered by the results presented 
here, because it is not clear how the genetic concept of fit- 
ness translates into demographic factors such as mosquito 
population size and longevity that determine the inten- 
sity of disease transmission. On the other hand, as note 
above, if synergist throws most of the selection pres- 
sure onto another niche then overall the rate of selection 
for resistance may increase. Consequently the public use 
of insecticide within the home (predominantly as wall 
sprays and/or bed nets) cannot be investigates isolated 
from other insecticide applications that mosquitoes may 
encounter during their lifetime. This suggest that the 
malaria community is correct in being alarmed at the 
often uncontrolled use of insecticides in applications such 
as agriculture [38,39]. 
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